home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / liboctave / CColVector.cc < prev    next >
C/C++ Source or Header  |  1997-03-07  |  17KB  |  841 lines
  1. // ColumnVector manipulations.
  2. /*
  3.  
  4. Copyright (C) 1996 John W. Eaton
  5.  
  6. This file is part of Octave.
  7.  
  8. Octave is free software; you can redistribute it and/or modify it
  9. under the terms of the GNU General Public License as published by the
  10. Free Software Foundation; either version 2, or (at your option) any
  11. later version.
  12.  
  13. Octave is distributed in the hope that it will be useful, but WITHOUT
  14. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16. for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with Octave; see the file COPYING.  If not, write to the Free
  20. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  21.  
  22. */
  23.  
  24. #if defined (__GNUG__)
  25. #pragma implementation
  26. #endif
  27.  
  28. #ifdef HAVE_CONFIG_H
  29. #include <config.h>
  30. #endif
  31.  
  32. #include <iostream.h>
  33.  
  34. #include "f77-fcn.h"
  35. #include "lo-error.h"
  36. #include "mx-base.h"
  37. #include "mx-inlines.cc"
  38. #include "oct-cmplx.h"
  39.  
  40. // Fortran functions we call.
  41.  
  42. extern "C"
  43. {
  44.   int F77_FCN (zgemv, ZGEMV) (const char*, const int&, const int&,
  45.                   const Complex&, const Complex*,
  46.                   const int&, const Complex*, const int&,
  47.                   const Complex&, Complex*, const int&,
  48.                   long);
  49. }
  50.  
  51. // Complex Column Vector class
  52.  
  53. ComplexColumnVector::ComplexColumnVector (const ColumnVector& a)
  54.    : MArray<Complex> (a.length ())
  55. {
  56.   for (int i = 0; i < length (); i++)
  57.     elem (i) = a.elem (i);
  58. }
  59.  
  60. bool
  61. ComplexColumnVector::operator == (const ComplexColumnVector& a) const
  62. {
  63.   int len = length ();
  64.   if (len != a.length ())
  65.     return 0;
  66.   return equal (data (), a.data (), len);
  67. }
  68.  
  69. bool
  70. ComplexColumnVector::operator != (const ComplexColumnVector& a) const
  71. {
  72.   return !(*this == a);
  73. }
  74.  
  75. // destructive insert/delete/reorder operations
  76.  
  77. ComplexColumnVector&
  78. ComplexColumnVector::insert (const ColumnVector& a, int r)
  79. {
  80.   int a_len = a.length ();
  81.   if (r < 0 || r + a_len > length ())
  82.     {
  83.       (*current_liboctave_error_handler) ("range error for insert");
  84.       return *this;
  85.     }
  86.  
  87.   for (int i = 0; i < a_len; i++)
  88.     elem (r+i) = a.elem (i);
  89.  
  90.   return *this;
  91. }
  92.  
  93. ComplexColumnVector&
  94. ComplexColumnVector::insert (const ComplexColumnVector& a, int r)
  95. {
  96.   int a_len = a.length ();
  97.   if (r < 0 || r + a_len > length ())
  98.     {
  99.       (*current_liboctave_error_handler) ("range error for insert");
  100.       return *this;
  101.     }
  102.  
  103.   for (int i = 0; i < a_len; i++)
  104.     elem (r+i) = a.elem (i);
  105.  
  106.   return *this;
  107. }
  108.  
  109. ComplexColumnVector&
  110. ComplexColumnVector::fill (double val)
  111. {
  112.   int len = length ();
  113.   if (len > 0)
  114.     for (int i = 0; i < len; i++)
  115.       elem (i) = val;
  116.   return *this;
  117. }
  118.  
  119. ComplexColumnVector&
  120. ComplexColumnVector::fill (const Complex& val)
  121. {
  122.   int len = length ();
  123.   if (len > 0)
  124.     for (int i = 0; i < len; i++)
  125.       elem (i) = val;
  126.   return *this;
  127. }
  128.  
  129. ComplexColumnVector&
  130. ComplexColumnVector::fill (double val, int r1, int r2)
  131. {
  132.   int len = length ();
  133.   if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
  134.     {
  135.       (*current_liboctave_error_handler) ("range error for fill");
  136.       return *this;
  137.     }
  138.  
  139.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  140.  
  141.   for (int i = r1; i <= r2; i++)
  142.     elem (i) = val;
  143.  
  144.   return *this;
  145. }
  146.  
  147. ComplexColumnVector&
  148. ComplexColumnVector::fill (const Complex& val, int r1, int r2)
  149. {
  150.   int len = length ();
  151.   if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
  152.     {
  153.       (*current_liboctave_error_handler) ("range error for fill");
  154.       return *this;
  155.     }
  156.  
  157.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  158.  
  159.   for (int i = r1; i <= r2; i++)
  160.     elem (i) = val;
  161.  
  162.   return *this;
  163. }
  164.  
  165. ComplexColumnVector
  166. ComplexColumnVector::stack (const ColumnVector& a) const
  167. {
  168.   int len = length ();
  169.   int nr_insert = len;
  170.   ComplexColumnVector retval (len + a.length ());
  171.   retval.insert (*this, 0);
  172.   retval.insert (a, nr_insert);
  173.   return retval;
  174. }
  175.  
  176. ComplexColumnVector
  177. ComplexColumnVector::stack (const ComplexColumnVector& a) const
  178. {
  179.   int len = length ();
  180.   int nr_insert = len;
  181.   ComplexColumnVector retval (len + a.length ());
  182.   retval.insert (*this, 0);
  183.   retval.insert (a, nr_insert);
  184.   return retval;
  185. }
  186.  
  187. ComplexRowVector
  188. ComplexColumnVector::hermitian (void) const
  189. {
  190.   int len = length ();
  191.   return ComplexRowVector (conj_dup (data (), len), len);
  192. }
  193.  
  194. ComplexRowVector
  195. ComplexColumnVector::transpose (void) const
  196. {
  197.   return ComplexRowVector (*this);
  198. }
  199.  
  200. ComplexColumnVector
  201. conj (const ComplexColumnVector& a)
  202. {
  203.   int a_len = a.length ();
  204.   ComplexColumnVector retval;
  205.   if (a_len > 0)
  206.     retval = ComplexColumnVector (conj_dup (a.data (), a_len), a_len);
  207.   return retval;
  208. }
  209.  
  210. // resize is the destructive equivalent for this one
  211.  
  212. ComplexColumnVector
  213. ComplexColumnVector::extract (int r1, int r2) const
  214. {
  215.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  216.  
  217.   int new_r = r2 - r1 + 1;
  218.  
  219.   ComplexColumnVector result (new_r);
  220.  
  221.   for (int i = 0; i < new_r; i++)
  222.     result.elem (i) = elem (r1+i);
  223.  
  224.   return result;
  225. }
  226.  
  227. // column vector by column vector -> column vector operations
  228.  
  229. ComplexColumnVector&
  230. ComplexColumnVector::operator += (const ColumnVector& a)
  231. {
  232.   int len = length ();
  233.  
  234.   int a_len = a.length ();
  235.  
  236.   if (len != a_len)
  237.     {
  238.       gripe_nonconformant ("operator +=", len, a_len);
  239.       return *this;
  240.     }
  241.  
  242.   if (len == 0)
  243.     return *this;
  244.  
  245.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  246.  
  247.   add2 (d, a.data (), len);
  248.   return *this;
  249. }
  250.  
  251. ComplexColumnVector&
  252. ComplexColumnVector::operator -= (const ColumnVector& a)
  253. {
  254.   int len = length ();
  255.  
  256.   int a_len = a.length ();
  257.  
  258.   if (len != a_len)
  259.     {
  260.       gripe_nonconformant ("operator -=", len, a_len);
  261.       return *this;
  262.     }
  263.  
  264.   if (len == 0)
  265.     return *this;
  266.  
  267.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  268.  
  269.   subtract2 (d, a.data (), len);
  270.   return *this;
  271. }
  272.  
  273. ComplexColumnVector&
  274. ComplexColumnVector::operator += (const ComplexColumnVector& a)
  275. {
  276.   int len = length ();
  277.  
  278.   int a_len = a.length ();
  279.  
  280.   if (len != a_len)
  281.     {
  282.       gripe_nonconformant ("operator +=", len, a_len);
  283.       return *this;
  284.     }
  285.  
  286.   if (len == 0)
  287.     return *this;
  288.  
  289.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  290.  
  291.   add2 (d, a.data (), len);
  292.   return *this;
  293. }
  294.  
  295. ComplexColumnVector&
  296. ComplexColumnVector::operator -= (const ComplexColumnVector& a)
  297. {
  298.   int len = length ();
  299.  
  300.   int a_len = a.length ();
  301.  
  302.   if (len != a_len)
  303.     {
  304.       gripe_nonconformant ("operator -=", len, a_len);
  305.       return *this;
  306.     }
  307.  
  308.   if (len == 0)
  309.     return *this;
  310.  
  311.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  312.  
  313.   subtract2 (d, a.data (), len);
  314.   return *this;
  315. }
  316.  
  317. // column vector by scalar -> column vector operations
  318.  
  319. ComplexColumnVector
  320. operator + (const ComplexColumnVector& v, double s)
  321. {
  322.   int len = v.length ();
  323.   return ComplexColumnVector (add (v.data (), len, s), len);
  324. }
  325.  
  326. ComplexColumnVector
  327. operator - (const ComplexColumnVector& v, double s)
  328. {
  329.   int len = v.length ();
  330.   return ComplexColumnVector (subtract (v.data (), len, s), len);
  331. }
  332.  
  333. ComplexColumnVector
  334. operator * (const ComplexColumnVector& v, double s)
  335. {
  336.   int len = v.length ();
  337.   return ComplexColumnVector (multiply (v.data (), len, s), len);
  338. }
  339.  
  340. ComplexColumnVector
  341. operator / (const ComplexColumnVector& v, double s)
  342. {
  343.   int len = v.length ();
  344.   return ComplexColumnVector (divide (v.data (), len, s), len);
  345. }
  346.  
  347. ComplexColumnVector
  348. operator + (const ColumnVector& a, const Complex& s)
  349. {
  350.   int len = a.length ();
  351.   return ComplexColumnVector (add (a.data (), len, s), len);
  352. }
  353.  
  354. ComplexColumnVector
  355. operator - (const ColumnVector& a, const Complex& s)
  356. {
  357.   int len = a.length ();
  358.   return ComplexColumnVector (subtract (a.data (), len, s), len);
  359. }
  360.  
  361. ComplexColumnVector
  362. operator * (const ColumnVector& a, const Complex& s)
  363. {
  364.   int len = a.length ();
  365.   return ComplexColumnVector (multiply (a.data (), len, s), len);
  366. }
  367.  
  368. ComplexColumnVector
  369. operator / (const ColumnVector& a, const Complex& s)
  370. {
  371.   int len = a.length ();
  372.   return ComplexColumnVector (divide (a.data (), len, s), len);
  373. }
  374.  
  375. // scalar by column vector -> column vector operations
  376.  
  377. ComplexColumnVector
  378. operator + (double s, const ComplexColumnVector& a)
  379. {
  380.   int a_len = a.length ();
  381.   return ComplexColumnVector (add (a.data (), a_len, s), a_len);
  382. }
  383.  
  384. ComplexColumnVector
  385. operator - (double s, const ComplexColumnVector& a)
  386. {
  387.   int a_len = a.length ();
  388.   return ComplexColumnVector (subtract (s, a.data (), a_len), a_len);
  389. }
  390.  
  391. ComplexColumnVector
  392. operator * (double s, const ComplexColumnVector& a)
  393. {
  394.   int a_len = a.length ();
  395.   return ComplexColumnVector (multiply (a.data (), a_len, s), a_len);
  396. }
  397.  
  398. ComplexColumnVector
  399. operator / (double s, const ComplexColumnVector& a)
  400. {
  401.   int a_len = a.length ();
  402.   return ComplexColumnVector (divide (s, a.data (), a_len), a_len);
  403. }
  404.  
  405. ComplexColumnVector
  406. operator + (const Complex& s, const ColumnVector& a)
  407. {
  408.   int a_len = a.length ();
  409.   return ComplexColumnVector (add (a.data (), a_len, s), a_len);
  410. }
  411.  
  412. ComplexColumnVector
  413. operator - (const Complex& s, const ColumnVector& a)
  414. {
  415.   int a_len = a.length ();
  416.   return ComplexColumnVector (subtract (s, a.data (), a_len), a_len);
  417. }
  418.  
  419. ComplexColumnVector
  420. operator * (const Complex& s, const ColumnVector& a)
  421. {
  422.   int a_len = a.length ();
  423.   return ComplexColumnVector (multiply (a.data (), a_len, s), a_len);
  424. }
  425.  
  426. ComplexColumnVector
  427. operator / (const Complex& s, const ColumnVector& a)
  428. {
  429.   int a_len = a.length ();
  430.   return ComplexColumnVector (divide (s, a.data (), a_len), a_len);
  431. }
  432.  
  433. // matrix by column vector -> column vector operations
  434.  
  435. ComplexColumnVector
  436. operator * (const ComplexMatrix& m, const ColumnVector& a)
  437. {
  438.   ComplexColumnVector tmp (a);
  439.   return m * tmp;
  440. }
  441.  
  442. ComplexColumnVector
  443. operator * (const ComplexMatrix& m, const ComplexColumnVector& a)
  444. {
  445.   ComplexColumnVector retval;
  446.  
  447.   int nr = m.rows ();
  448.   int nc = m.cols ();
  449.  
  450.   int a_len = a.length ();
  451.  
  452.   if (nc != a_len)
  453.     gripe_nonconformant ("operator *", nr, nc, a_len, 1);
  454.   else
  455.     {
  456.       if (nc == 0 || nr == 0)
  457.     retval.resize (nr, 0.0);
  458.       else
  459.     {
  460.       int ld = nr;
  461.  
  462.       retval.resize (nr);
  463.       Complex *y = retval.fortran_vec ();
  464.  
  465.       F77_XFCN (zgemv, ZGEMV, ("N", nr, nc, 1.0, m.data (), ld,
  466.                    a.data (), 1, 0.0, y, 1, 1L));
  467.  
  468.       if (f77_exception_encountered)
  469.         (*current_liboctave_error_handler)
  470.           ("unrecoverable error in zgemv");
  471.     }
  472.     }
  473.  
  474.   return retval;
  475. }
  476.  
  477. // column vector by column vector -> column vector operations
  478.  
  479. ComplexColumnVector
  480. operator + (const ComplexColumnVector& v, const ColumnVector& a)
  481. {
  482.   int len = v.length ();
  483.  
  484.   int a_len = a.length ();
  485.  
  486.   if (len != a_len)
  487.     {
  488.       gripe_nonconformant ("operator +", len, a_len);
  489.       return ComplexColumnVector ();
  490.     }
  491.  
  492.   if (len == 0)
  493.     return ComplexColumnVector (0);
  494.  
  495.   return ComplexColumnVector (add (v.data (), a.data (), len), len);
  496. }
  497.  
  498. ComplexColumnVector
  499. operator - (const ComplexColumnVector& v, const ColumnVector& a)
  500. {
  501.   int len = v.length ();
  502.  
  503.   int a_len = a.length ();
  504.  
  505.   if (len != a_len)
  506.     {
  507.       gripe_nonconformant ("operator -", len, a_len);
  508.       return ComplexColumnVector ();
  509.     }
  510.  
  511.   if (len == 0)
  512.     return ComplexColumnVector (0);
  513.  
  514.   return ComplexColumnVector (subtract (v.data (), a.data (), len), len);
  515. }
  516.  
  517. ComplexColumnVector
  518. operator + (const ColumnVector& v, const ComplexColumnVector& a)
  519. {
  520.   int len = v.length ();
  521.  
  522.   int a_len = a.length ();
  523.  
  524.   if (len != a_len)
  525.     {
  526.       gripe_nonconformant ("operator +", len, a_len);
  527.       return ComplexColumnVector ();
  528.     }
  529.  
  530.   if (len == 0)
  531.     return ComplexColumnVector (0);
  532.  
  533.   return ComplexColumnVector (add (v.data (), a.data (), len), len);
  534. }
  535.  
  536. ComplexColumnVector
  537. operator - (const ColumnVector& v, const ComplexColumnVector& a)
  538. {
  539.   int len = v.length ();
  540.  
  541.   int a_len = a.length ();
  542.  
  543.   if (len != a_len)
  544.     {
  545.       gripe_nonconformant ("operator -", len, a_len);
  546.       return ComplexColumnVector ();
  547.     }
  548.  
  549.   if (len == 0)
  550.     return ComplexColumnVector (0);
  551.  
  552.   return ComplexColumnVector (subtract (v.data (), a.data (), len), len);
  553. }
  554.  
  555. ComplexColumnVector
  556. product (const ComplexColumnVector& v, const ColumnVector& a)
  557. {
  558.   int len = v.length ();
  559.  
  560.   int a_len = a.length ();
  561.  
  562.   if (len != a_len)
  563.     {
  564.       gripe_nonconformant ("product", len, a_len);
  565.       return ComplexColumnVector ();
  566.     }
  567.  
  568.   if (len == 0)
  569.     return ComplexColumnVector (0);
  570.  
  571.   return ComplexColumnVector (multiply (v.data (), a.data (), len), len);
  572. }
  573.  
  574. ComplexColumnVector
  575. quotient (const ComplexColumnVector& v, const ColumnVector& a)
  576. {
  577.   int len = v.length ();
  578.  
  579.   int a_len = a.length ();
  580.  
  581.   if (len != a_len)
  582.     {
  583.       gripe_nonconformant ("quotient", len, a_len);
  584.       return ComplexColumnVector ();
  585.     }
  586.  
  587.   if (len == 0)
  588.     return ComplexColumnVector (0);
  589.  
  590.   return ComplexColumnVector (divide (v.data (), a.data (), len), len);
  591. }
  592.  
  593. ComplexColumnVector
  594. product (const ColumnVector& v, const ComplexColumnVector& a)
  595. {
  596.   int len = v.length ();
  597.  
  598.   int a_len = a.length ();
  599.  
  600.   if (len != a_len)
  601.     {
  602.       gripe_nonconformant ("product", len, a_len);
  603.       return ColumnVector ();
  604.     }
  605.  
  606.   if (len == 0)
  607.     return ComplexColumnVector (0);
  608.  
  609.   return ComplexColumnVector (multiply (v.data (), a.data (), len), len);
  610. }
  611.  
  612. ComplexColumnVector
  613. quotient (const ColumnVector& v, const ComplexColumnVector& a)
  614. {
  615.   int len = v.length ();
  616.  
  617.   int a_len = a.length ();
  618.  
  619.   if (len != a_len)
  620.     {
  621.       gripe_nonconformant ("quotient", len, a_len);
  622.       return ColumnVector ();
  623.     }
  624.  
  625.   if (len == 0)
  626.     return ComplexColumnVector (0);
  627.  
  628.   return ComplexColumnVector (divide (v.data (), a.data (), len), len);
  629. }
  630.  
  631. // matrix by column vector -> column vector operations
  632.  
  633. ComplexColumnVector
  634. operator * (const Matrix& m, const ComplexColumnVector& a)
  635. {
  636.   ComplexMatrix tmp (m);
  637.   return tmp * a;
  638. }
  639.  
  640. // diagonal matrix by column vector -> column vector operations
  641.  
  642. ComplexColumnVector
  643. operator * (const DiagMatrix& m, const ComplexColumnVector& a)
  644. {
  645.   int nr = m.rows ();
  646.   int nc = m.cols ();
  647.  
  648.   int a_len = a.length ();
  649.  
  650.   if (nc != a_len)
  651.     {
  652.       gripe_nonconformant ("operator *", nr, nc, a_len, 1);
  653.       return ColumnVector ();
  654.     }
  655.  
  656.   if (nc == 0 || nr == 0)
  657.     return ComplexColumnVector (0);
  658.  
  659.   ComplexColumnVector result (nr);
  660.  
  661.   for (int i = 0; i < a_len; i++)
  662.     result.elem (i) = a.elem (i) * m.elem (i, i);
  663.  
  664.   for (int i = a_len; i < nr; i++)
  665.     result.elem (i) = 0.0;
  666.  
  667.   return result;
  668. }
  669.  
  670. ComplexColumnVector
  671. operator * (const ComplexDiagMatrix& m, const ColumnVector& a)
  672. {
  673.   int nr = m.rows ();
  674.   int nc = m.cols ();
  675.  
  676.   int a_len = a.length ();
  677.  
  678.   if (nc != a_len)
  679.     {
  680.       gripe_nonconformant ("operator *", nr, nc, a_len, 1);
  681.       return ComplexColumnVector ();
  682.     }
  683.  
  684.   if (nc == 0 || nr == 0)
  685.     return ComplexColumnVector (0);
  686.  
  687.   ComplexColumnVector result (nr);
  688.  
  689.   for (int i = 0; i < a_len; i++)
  690.     result.elem (i) = a.elem (i) * m.elem (i, i);
  691.  
  692.   for (int i = a_len; i < nr; i++)
  693.     result.elem (i) = 0.0;
  694.  
  695.   return result;
  696. }
  697.  
  698. ComplexColumnVector
  699. operator * (const ComplexDiagMatrix& m, const ComplexColumnVector& a)
  700. {
  701.   int nr = m.rows ();
  702.   int nc = m.cols ();
  703.  
  704.   int a_len = a.length ();
  705.  
  706.   if (nc != a_len)
  707.     {
  708.       gripe_nonconformant ("operator *", nr, nc, a_len, 1);
  709.       return ComplexColumnVector ();
  710.     }
  711.  
  712.   if (nc == 0 || nr == 0)
  713.     return ComplexColumnVector (0);
  714.  
  715.   ComplexColumnVector result (nr);
  716.  
  717.   for (int i = 0; i < a_len; i++)
  718.     result.elem (i) = a.elem (i) * m.elem (i, i);
  719.  
  720.   for (int i = a_len; i < nr; i++)
  721.     result.elem (i) = 0.0;
  722.  
  723.   return result;
  724. }
  725.  
  726. // other operations
  727.  
  728. ComplexColumnVector
  729. ComplexColumnVector::map (c_c_Mapper f) const
  730. {
  731.   ComplexColumnVector b (*this);
  732.   return b.apply (f);
  733. }
  734.  
  735. ColumnVector
  736. ComplexColumnVector::map (d_c_Mapper f) const
  737. {
  738.   const Complex *d = data ();
  739.  
  740.   int len = length ();
  741.  
  742.   ColumnVector retval (len);
  743.  
  744.   double *r = retval.fortran_vec ();
  745.  
  746.   for (int i = 0; i < len; i++)
  747.     r[i] = f (d[i]);
  748.  
  749.   return retval;
  750. }
  751.  
  752. ComplexColumnVector&
  753. ComplexColumnVector::apply (c_c_Mapper f)
  754. {
  755.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  756.  
  757.   for (int i = 0; i < length (); i++)
  758.     d[i] = f (d[i]);
  759.  
  760.   return *this;
  761. }
  762.  
  763. Complex
  764. ComplexColumnVector::min (void) const
  765. {
  766.   int len = length ();
  767.   if (len == 0)
  768.     return 0.0;
  769.  
  770.   Complex res = elem (0);
  771.   double absres = abs (res);
  772.  
  773.   for (int i = 1; i < len; i++)
  774.     if (abs (elem (i)) < absres)
  775.       {
  776.     res = elem (i);
  777.     absres = abs (res);
  778.       }
  779.  
  780.   return res;
  781. }
  782.  
  783. Complex
  784. ComplexColumnVector::max (void) const
  785. {
  786.   int len = length ();
  787.   if (len == 0)
  788.     return 0.0;
  789.  
  790.   Complex res = elem (0);
  791.   double absres = abs (res);
  792.  
  793.   for (int i = 1; i < len; i++)
  794.     if (abs (elem (i)) > absres)
  795.       {
  796.     res = elem (i);
  797.     absres = abs (res);
  798.       }
  799.  
  800.   return res;
  801. }
  802.  
  803. // i/o
  804.  
  805. ostream&
  806. operator << (ostream& os, const ComplexColumnVector& a)
  807. {
  808. //  int field_width = os.precision () + 7;
  809.   for (int i = 0; i < a.length (); i++)
  810.     os << /* setw (field_width) << */ a.elem (i) << "\n";
  811.   return os;
  812. }
  813.  
  814. istream&
  815. operator >> (istream& is, ComplexColumnVector& a)
  816. {
  817.   int len = a.length();
  818.  
  819.   if (len < 1)
  820.     is.clear (ios::badbit);
  821.   else
  822.     {
  823.       double tmp;
  824.       for (int i = 0; i < len; i++)
  825.         {
  826.           is >> tmp;
  827.           if (is)
  828.             a.elem (i) = tmp;
  829.           else
  830.             break;
  831.         }
  832.     }
  833.   return is;
  834. }
  835.  
  836. /*
  837. ;;; Local Variables: ***
  838. ;;; mode: C++ ***
  839. ;;; End: ***
  840. */
  841.